#include "dat3d.h"
#include <fstream>
#include <cstdlib>
#include "string.h"
//#include <exception>


#ifndef n_MIN
	#define n_MIN (double)0.01 ///Mpc
#endif
#define l_MIN log10(n_MIN)

#ifndef n_MAX 
	#define n_MAX (double)50.0
#endif
#define l_MAX log10(n_MAX)


/*#ifndef _BINES
	#define _BINES 20
#endif*/

#ifndef n_PIZZA
	#define n_PIZZA (double) 10
#endif

#ifndef _max_Z
	#define _max_Z 0.10
#endif
#ifndef _min_Z
	#define _min_Z 0.02
#endif
#define N_JACK 8
// los pongo afuera porque lleve leemos afuera.... hay que ver si conviene o los pasamos entre las funciones (hacer una union?)
int rap=-1;
Healpix_Base grid(8,RING);// grid
//vector<int> fin[grid.Npix()];	// Indices de las galalaxias en funcion del bin de healpix
void help(void);
void ran_ini(void);

int main(int argc,char * argv[])
{
	vector<pointing> p_ang,p_ang_trac; //pointing_angle
	vector<double> Z,Z_trac; //Zs
	vector<double> marc,marc_trac,mtmp; // marquita
	vector<long int> hpx_idx;
	
	vector<int> *fin;	// Indices de las galalaxias en funcion del bin de healpix
	vector<int> *fin_trac;	// Indices de las galalaxias en funcion del bin de healpix
	fin=(vector<int> * )malloc(sizeof(vector<int>)*grid.Npix());
	fin_trac=(vector<int> * )malloc(sizeof(vector<int>)*grid.Npix());

	double sig,pi;
	int cual,cualp,cualrs,i,j,k,entero,m,n,ii;
	double d_tmp,d_tmp_bis;

	string ss = "salida.dat";
	string n_dat = "entrada";
	string n_dat_trac = "entrada";
	string mensajes = "CALCULANDO  ";
	
	uni_type U1;
	cout << "-------------------VERSION "<< VERSION <<" -------------------------\n";
	cout << "*****************************************************\n";
	cout <<" Npix : "<< grid.Npix() << "\n";fflush(stdout);	
	for(i=1;i<argc;i++){
               	if(!strcmp(argv[i],"-s")){ 
			if(!ss.compare("salida.dat")) ss.assign(argv[++i]);
			cout << "OK salida Marcos: " << ss << endl;
                }else if(!strcmp(argv[i],"-Nd")){
                        if(!n_dat.compare("entrada")) n_dat.assign(argv[++i]);
			cout << "OK! la entrada 1: " << n_dat << endl;
			if(i+1<argc)
			if(argv[i+1][0]!='-'){
				n_dat_trac.assign(argv[++i]);
				cout << "OK! la entrada 2: " << n_dat_trac << endl;
			}
                }else if(!strcmp(argv[i],"-Xrap")){
                        if(rap==-1) rap=atoi(argv[++i]);
                        printf("OK! BOOST: %i \n", rap);
                }else if(!strcmp(argv[i],"-help")){
			help();
		}else{ 
			cout << "\nparametro " << argv[i] << " no funca. "<< endl;
		       	cout << "ie NO SIRVE PA NADA SACALO DE AHI!!!\n";
		       	cout << "intente :\n " <<argv[0] <<" -help" << endl;
			exit(2);
		}
	};
	U1=universo(0.3,0.7,_max_Z);

	cout << "leemos centros" << endl;
	leemos(n_dat,Z,p_ang,marc,fin);
	cout << "leemos tracers" << endl;
	leemos(n_dat_trac,Z_trac,p_ang_trac,marc_trac,fin_trac);
	
	double m12=0.0;  // esl m(r) medio que uso en la fromula
	double a[3],b[3];   // vectores temporales para el calculo de distancia
	double lola,lola2;  // angulo de poscicion  y 
#ifdef _BINES
	int bines=_BINES;   
#else
	int bines=(int)floor(sqrt(Z.size()));   
#endif
	long double count[bines];
	long double count_NUM[bines];
	long double count_err[N_JACK][bines];
	long double count_err1[N_JACK][bines];
	long double count_err2[N_JACK][bines];
	long double count_err_med[bines];
	long int    count_count[N_JACK][bines];
	long double boroto[N_JACK][bines][bines];
	long double boroto_rand[N_JACK][bines][bines];
	long double boroto_u[N_JACK][bines][bines];
	long double boroto_u_rand[N_JACK][bines][bines];
	long double boroto_DD[bines][bines];
	long double boroto_count[bines][bines];
	long double boroto_DR[bines][bines];
	long double boroto_RR[bines][bines];
	long int kk,jj;
	double memedia=0, test_0,test_1;
	vector <int> interesantes;   //indices de las glx que estan en X bin
	vector <int> menosinteresantes;   //para hacer un ring?
	vector <int>::iterator it,itt;
	vector <int> j_rand,i_rand,k_rand;
	vector <int> i_jack;
	vector <double> xyz[3];
	vector <double> xyz_trac[3];

	int deb,debb,debbb;  ///mas debuggeo
	deb=debb=debbb=0;
	

	for(i=0;i<bines;i++) {count[i]=0.; count_NUM[i]=0.;}
	for(i=0;i<bines;i++) count_err_med[i]=0.;
	for(i=0;i<bines;i++) for(jj=0;jj<N_JACK;jj++) count_err[jj][i]=0.;
	for(i=0;i<bines;i++) for(jj=0;jj<N_JACK;jj++) count_err1[jj][i]=0.;
	for(i=0;i<bines;i++) for(jj=0;jj<N_JACK;jj++) count_err2[jj][i]=0.;
	for(i=0;i<bines;i++) for(jj=0;jj<N_JACK;jj++) count_count[jj][i]=0;
	for(i=0;i<bines;i++) for(k=0;k<bines;k++) for(jj=0;jj<N_JACK;jj++){ 
		boroto[jj][i][k]=0.; boroto_rand[jj][i][k]=0.;
		boroto_u[jj][i][k]=0.; boroto_u_rand[jj][i][k]=0.;}
	for(i=0;i<bines;i++) for(k=0;k<bines;k++) { boroto_DD[i][k]=0.; boroto_DR[i][k]=0.;boroto_RR[i][k]=0.;boroto_count[i][k]=0.;}
	
	cout << "Preparando datos auxiliares....";
	ran_ini();
	test_0=0;
	test_1=4*M_PI;
	for(i=0;i<Z.size();i++)
	{
		a[0]=sin(p_ang[i].theta)*cos(p_ang[i].phi);
		a[1]=sin(p_ang[i].theta)*sin(p_ang[i].phi);
		a[2]=cos(p_ang[i].theta);
		xyz[0].push_back(a[0]);
		xyz[1].push_back(a[1]);
		xyz[2].push_back(a[2]);
		test_0=p_ang[i].theta > test_0 ? p_ang[i].theta : test_0;
		test_1=p_ang[i].theta < test_1 ? p_ang[i].theta : test_1;

	//	i_jack.push_back((int)floor((double)N_JACK*(double)i/(double)Z.size()));
	}
	for(i=0;i<Z.size();i++)
	{	
		i_jack.push_back((int)floor((double) (p_ang[i].theta-test_1)/(test_0-test_1)*N_JACK));	
	}

	cout << endl << " Centros listos " << Z.size() << endl;
	
	for(i=0;i<Z_trac.size();i++)
	{
		a[0]=sin(p_ang_trac[i].theta)*cos(p_ang_trac[i].phi);
		a[1]=sin(p_ang_trac[i].theta)*sin(p_ang_trac[i].phi);
		a[2]=cos(p_ang_trac[i].theta);
		xyz_trac[0].push_back(a[0]);
		xyz_trac[1].push_back(a[1]);
		xyz_trac[2].push_back(a[2]);
		j_rand.push_back(i);
		i_rand.push_back(i);
		k_rand.push_back(i);
	}
	
	for(i=0;i<Z_trac.size();i++){
//		jj=j_rand[i];
		kk=(int)floor(drand48()*(double)Z_trac.size());
		j_rand[kk]=i;
		j_rand[i]=kk;
		kk=(int)floor(drand48()*(double)Z_trac.size());
		i_rand[kk]=i;
		i_rand[i]=kk;
		kk=(int)floor(drand48()*(double)Z_trac.size());
		k_rand[kk]=i;
		k_rand[i]=kk;
	}
	
	cout << " Tracer  listos "<<Z_trac.size() << endl;

	for(i=0;i<Z.size();i++)
	{
	    d_tmp=asin(n_MAX/(rr(U1,Z[i])));
	    grid.query_disc_inclusive(p_ang[i],d_tmp,interesantes);
	    d_tmp=asin(n_MIN/(rr(U1,Z[i])));
	    grid.query_disc(p_ang[i],d_tmp,menosinteresantes);
	   
	      for(kk=0;kk<menosinteresantes.size();kk++){
	      	it=remove(interesantes.begin(),interesantes.end(),menosinteresantes[kk]);
	        for(itt=interesantes.end();itt!=it;itt--)
		  interesantes.pop_back();
	      }

	      a[0]=xyz[0][i];
	      a[1]=xyz[1][i];
	      a[2]=xyz[2][i];
	      
	      for(k=0;k<interesantes.size();k++){ // boocle sobre los bines de Healpix que nos interesan
		  for(kk=0;kk<fin_trac[interesantes[k]].size();kk++) 
		  // boocle sobre las galaxias dentro del Bin de interes
		  {
		     j=fin_trac[interesantes[k]][kk];
		//DATA - DATA
		     if(abs(rr(U1,Z[i])-rr(U1,Z_trac[j]))<abs(rr(U1,Z[i])-(n_PIZZA)))
		     { 	// intento hacerlo mas rapico cortando los Z que no gustan 
			b[0]=xyz_trac[0][j];
			b[1]=xyz_trac[1][j];
			b[2]=xyz_trac[2][j];
			
			 d_tmp_bis=acos((a[0]*b[0]+a[1]*b[1]+a[2]*b[2]));
			 
			 sig=log10((double)(rr(U1,Z_trac[j])+rr(U1,Z[i]))*sin(d_tmp_bis/2));
			 pi= log10((double)fabs(rr(U1,Z_trac[j])-rr(U1,Z[i]))*cos(d_tmp_bis/2));

//			lola=log10(dist(a,b));
//			lola=log10(acos((a[0]*b[0]+a[1]*b[1]+a[2]*b[2]))*rr(U1,Z[i]));
                        
			if((pi>l_MIN && pi<=l_MAX) &&(sig>l_MIN && sig<=l_MAX)){
				cual= (int)floor((sig-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				cualp= (int)floor((pi-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				boroto[i_jack[i]][cual][cualp]++;			
			}
		     }
		// DATA - RANDOM
		     jj=j_rand[j];
//		     cout << jj << Z_trac.size();
		     if(abs(rr(U1,Z[i])-rr(U1,Z_trac[jj]))<abs(rr(U1,Z[i])-(n_PIZZA)))
		     { 	// intento hacerlo mas rapico cortando los Z que no gustan 
			 
			 sig=log10((double)(rr(U1,Z_trac[jj])+rr(U1,Z[i]))*sin(d_tmp_bis/2));
			 pi= log10((double)fabs(rr(U1,Z_trac[jj])-rr(U1,Z[i]))*cos(d_tmp_bis/2));

			if((pi>l_MIN && pi<=l_MAX) && (sig>l_MIN && sig<=l_MAX)){
				cual = (int)floor((sig-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				cualp= (int)floor((pi-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				boroto_rand[i_jack[i]][cual][cualp]++;		
			}
		     }
		     jj=k_rand[j];
//		     cout << jj << Z_trac.size();
		     if(abs(rr(U1,Z[i])-rr(U1,Z_trac[jj]))<abs(rr(U1,Z[i])-(n_PIZZA)))
		     { 	// intento hacerlo mas rapico cortando los Z que no gustan 
			 
			 sig=log10((double)(rr(U1,Z_trac[jj])+rr(U1,Z[i]))*sin(d_tmp_bis/2));
			 pi= log10((double)fabs(rr(U1,Z_trac[jj])-rr(U1,Z[i]))*cos(d_tmp_bis/2));

			if((pi>l_MIN && pi<=l_MAX) && (sig>l_MIN && sig<=l_MAX)){
				cual = (int)floor((sig-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				cualp= (int)floor((pi-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				boroto_rand[i_jack[i]][cual][cualp]++;		
			}
		     }
		     jj=i_rand[j];
//		     cout << jj << Z_trac.size();
		     if(abs(rr(U1,Z[i])-rr(U1,Z_trac[jj]))<abs(rr(U1,Z[i])-(n_PIZZA)))
		     { 	// intento hacerlo mas rapico cortando los Z que no gustan 
			 
			 sig=log10((double)(rr(U1,Z_trac[jj])+rr(U1,Z[i]))*sin(d_tmp_bis/2));
			 pi= log10((double)fabs(rr(U1,Z_trac[jj])-rr(U1,Z[i]))*cos(d_tmp_bis/2));

			if((pi>l_MIN && pi<=l_MAX) && (sig>l_MIN && sig<=l_MAX)){
				cual = (int)floor((sig-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				cualp= (int)floor((pi-l_MIN)/(l_MAX-l_MIN)*bines);	// busco el bin
				boroto_rand[i_jack[i]][cual][cualp]++;		
			}
		     }
			//fin conteo
		  }
	      }  // fin interesantes
	    if(rap<0)
	    graf2(&mensajes[0],i,Z.size(),45);
	}
	
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
/*	for(i=0;i<bines;i++)
		for(jj=0;jj<N_JACK;jj++){
		count[i]+=count_err[jj][i];
		count_NUM[i]+=count_count[jj][i];
	}
	for(i=0;i<bines;i++)
		count[i]=(double)count[i]/(median*median_trac*(double)count_NUM[i]);
	for(i=0;i<bines;i++)
		for(kk=0;kk<N_JACK;kk++)
		for(jj=0;jj<N_JACK;jj++){
			count_err1[jj][i]+=count_err[kk][i];//(median*median_trac*count_count[kk][i])/(N_JACK-1);
			count_err2[jj][i]+=count_count[kk][i];
		}
	for(i=0;i<bines;i++)
		for(jj=0;jj<N_JACK;jj++){
			count_err1[jj][i]-=count_err[jj][i];
			count_err2[jj][i]-=count_count[jj][i];
	}
	for(i=0;i<bines;i++)
		for(jj=0;jj<N_JACK;jj++)
			count_err[jj][i]=count_err1[jj][i]/(median*median_trac*count_err2[jj][i]);
	for(i=0;i<bines;i++)
		for(jj=0;jj<N_JACK;jj++){
			count_err_med[i]+=pow((double) (count[i]-count_err[jj][i]),2.);
		}	
	for(i=0;i<bines;i++)
		count_err_med[i]=sqrt(N_JACK*count_err_med[i]/(N_JACK-1));
*/
	cout << endl << "Starting statistics ... " << endl;
	
	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++)
	      for(k=0;k<N_JACK;k++){
	      	 boroto_DD[i][j]+=boroto[k][i][j];
	      	 boroto_DR[i][j]+=boroto_rand[k][i][j];
		}

	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++){
	   	if(boroto_DR[i][j]!=0){
	   	boroto_DD[i][j]=3.*boroto_DD[i][j]/boroto_DR[i][j]-1.; // aca esta la cantidad de veces que hay mas de randoms que de data
		boroto_DR[i][j]=0.;
		}
		else	
		boroto_DD[i][j]=0.;
	   }
	// este fue el verdadero
	// ahora los Jacks
	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++)
	      for(jj=0;jj<N_JACK;jj++)
	      for(k=0;k<N_JACK;k++){
	      	if(jj != k){
	      	      boroto_u[jj][i][j]=+boroto[k][i][j];
	      	 boroto_u_rand[jj][i][j]=+boroto_rand[k][i][j];}
		}
/*	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++)
	      for(k=0;k<N_JACK;k++){
	      	      boroto_u[k][i][j]-=boroto[k][i][j];
	      	 boroto_u_rand[k][i][j]-=boroto_rand[k][i][j];
		}*/
	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++){
	      for(k=0;k<N_JACK;k++){
	      		if(boroto_u_rand[k][i][j]!=0){
	      		boroto_u[k][i][j]=3.*boroto_u[k][i][j]/boroto_u_rand[k][i][j]-1.;//aca de nuevo
			boroto_count[i][j]++;
			}
			else
	      		boroto_u[k][i][j]=0.;//boroto_DD[i][j];

			boroto_RR[i][j]+=boroto_u[k][i][j];
		}
	      
	      }

	for(i=0;i<bines;i++)
	   for(j=0;j<bines;j++)
//	      for(k=0;k<N_JACK;k++)
	      	boroto_RR[i][j]/=boroto_count[i][j];
	      		
	for(i=0;i<bines;i++)
	for(j=0;j<bines;j++)
		for(jj=0;jj<N_JACK;jj++){
			if(boroto_u[jj][i][j] != 0.)
			boroto_DR[i][j]+=pow((double)(boroto_RR[i][j]-boroto_u[jj][i][j]),2.);
			else
			boroto_DR[i][j]+=0.;
		}	

	for(i=0;i<bines;i++)
	for(j=0;j<bines;j++)
		if(boroto_count[i][j] != 0)
		//boroto_DR[i][j]=sqrt(boroto_count[i][j]*boroto_DR[i][j]/(boroto_count[i][j]-1));
		boroto_DR[i][j]=sqrt(boroto_DR[i][j]/(boroto_count[i][j]-1));
		else
		boroto_DR[i][j]=0.;
//		boroto_DR[i][j]=sqrt(N_JACK*boroto_DR[i][j]/(N_JACK-1));
//		boroto_DR[i][j]=sqrt(boroto_DR[i][j]/(N_JACK-1));
	
	// terminamos con esto
	//
	printf("\nEscribimos\n");fflush(stdout);
	cout << "BINES: " << bines << " Ver.: "<< VERSION <<endl;	
	fflush(stdout);
	
	ofstream salida(&ss[0],ios::out);if(!salida){printf("error en %s para salida\n",&ss[0]); exit (7);};
	salida << "BINES: " << bines << " Ver.: "<< VERSION <<endl;	
	for(i=0;i<bines;i++)
	{
		d_tmp_bis=((2*i+1)*(l_MAX-l_MIN)/bines+2*l_MIN)/2;

		salida 	<< pow(10,d_tmp_bis) << endl; 
	}
	for(i=0;i<bines;i++)//sigma
	for(j=0;j<bines;j++)//PI
	{
		salida << boroto_DD[i][j] << " " <<boroto_DR[i][j] << " " << boroto_RR[i][j] << " " << boroto_count[i][j] << endl;
	}
	
	salida.close();salida.clear();
	cout << "lista escritura\n";
	printf("\nsalimos de 10\n");	
	return 0;
}

int leemos(string n_dat,vector<double> &Z_out,vector<pointing> &point,vector<double> &marc_out,vector<int> *fin)
{
	ifstream ent;
	int i;
	double dum_Z,dum_tmp;
	pointing dum_point;
	ent.open(&n_dat[0]);if(!ent){printf("error en %s\n",&n_dat[0]); exit (7);};
	for(i=0;!ent.eof();){
		ent >> dum_point.phi >> dum_point.theta >> dum_tmp >>  dum_Z;
		dum_point.phi=M_PI*dum_point.phi/180;
		dum_point.theta=M_PI*(90-dum_point.theta)/180;
		if(dum_Z>_min_Z && dum_Z<_max_Z ){
			point.push_back(dum_point);
			point[i].theta=dum_point.theta;
			point[i].phi=dum_point.phi;
			Z_out.push_back(dum_Z);
			marc_out.push_back(exp(dum_tmp));
			if(point[i].theta<0 || point[i].theta>M_PI){
				cerr << i <<" prolemas en angulos thita "<< endl; exit (9);}
			if(point[i].phi<0 || point[i].phi>2*M_PI){
				cerr << "prolemas en angulos phi "<< endl; exit (9);}
			fin[grid.ang2pix(point[i])].push_back(i);
			if(rap<0)
			graf2(&n_dat[0],i,100,30);
			i++;
		}
	}
	ent.close();ent.clear();
	cout << "\n";
	dum_Z=0;
	for (i=0;i<grid.Npix();i++){
		dum_Z=dum_Z > fin[i].size() ? dum_Z : fin[i].size();
		dum_tmp+=fin[i].size()/grid.Npix();
	}
	cout << "Maximo bin de healpix con: " << dum_Z << endl;
	return 0;
}

void help(void)
{
	cout << "****************************************************************" << endl;
	cout << "ok, " << endl;
	cout << "-Nd entra el acrhivo ascii" << endl;
	cout << "-s archivo de salida " << endl;
	exit(0);
};

uni_type universo(float omega_m,float omega_l,float Zmax,int iter)
{
  //lo que quiero hacer r(z)=c/Ho intergral 1/E(Z)
  int i;
  float b;
  uni_type resul;
  resul.binz=iter;
  resul.Z=new float[iter];
  resul.R=new float[iter];
  resul.Z[0]=0.;
  resul.R[0]=0.;
  cout << " Generando parametros de Universo" << endl;
  for(i = 1 ;i < iter; i++ ){
  	b=Zmax*float(i)/float(iter-1);
	resul.Z[i]=b;
	resul.R[i]=resul.R[i-1]+Z2R*(b-resul.Z[i-1])/sqrt(omega_m*pow(1+b,3)+omega_l);
//	cout << resul.Z[i] << "  " << resul.R[i] << " vs " << Z2R*b<< endl;
//	cout << b << " " << i << " " << Zmax << endl;
  } 
  cout << " Universo finalizado" << endl;
  return resul;
}

double rr(uni_type uni,float z)
{ 
  int i;
  if(uni.Z[uni.binz-1]<z){
        cout << "error en Z. Fuera del rango del universo" << endl;
  	exit(9090);
  }
  i=uni.binz-1;
  i=(int)floor(z/uni.Z[i]*i);
  return double(z*(uni.R[i+1]-uni.R[i])/(uni.Z[i+1]-uni.Z[i]));
};

void ran_ini(void)
{ ifstream ent; 
  long int seedval;

  ent.open("/dev/urandom",ios::binary);
  if(!ent){printf("error en INICIALIZAR RANDOMS"); exit (8);};
        ent.read((char *) & seedval,sizeof(seedval));
        ent.close();ent.clear();
        srand48(seedval);
//        cout << "SEED VAL: " << seedval << endl;
}
